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A unified FETI-DP approach for incompressible 

Stokes equations 

Xuemin Tu* Jing Lit 

Abstract 

A unified framework of FETI-DP algorithms is proposed for solving 
the system of linear equations arising from the mixed finite element ap- 
proximation of incompressible Stokes equations. A distinctive feature of 
this framework is that it allows using both continuous and discontinu- 
ous pressures in the algorithm, while previous FETI-DP methods only 
apply to discontinuous pressures. A preconditioned conjugate gradient 
method is used in the algorithm with either a lumped or a Dirichlet pre- 
conditioner, and scalable convergence rates are proved. This framework 
is also used to describe several previously developed FETI-DP algorithms 
and greatly simplifies their analysis. Numerical experiments of solving a 
two-dimensional incompressible Stokes problem demonstrate the perfor- 
mances of the discussed FETI-DP algorithms represented under the same 
framework. 

1 Introduction 

The finite element tearing and interconnecting (FETI) methods were introduced 
by Far hat and Roux [9] [TUJ [TT] [12] for second order elliptic problems. In these 
algorithms, Lagrange multipliers are introduced to enforce the continuity of the 
solution across the subdomain interface. The original system of linear equations 
is reduced to a symmetric positive semi-definite system for the Lagrange mul- 
tipliers, which can be solved by a preconditioned conjugate gradient method. 
Both a lumped preconditioner |10j and a Dirichlet preconditioner |12) have been 
used in the FETI methods. Compared with the lumped preconditioner, the 
Dirichlet preconditioner is more effective in the reduction of iteration count, but 
it is also more expensive. Numerical experiments in |12] show that the lumped 
preconditioner is more efficient for second-order problems while the Dirichlet 
preconditioner is better for plate and shell problems. 

The dual-primal FETI (FETI-DP) method, introduced by Farhat et. al. 
[7J IB], represents a further development of the FETI methods. In a FETI-DP 
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algorithm, a few variables from each subdomain are selected as the coarse level 
primal variables which are shared by neighboring subdomains, while the conti- 
nuity of the other subdomain interface variables is enforced by using Lagrange 
multipliers. The reduced system for the Lagrange multipliers becomes sym- 
metric positive definite and is solved by a preconditioned conjugate gradient 
method using either the lumped or the Dirichlet preconditioner. With an ap- 
propriate choice of coarse level variables, the condition number bound of the 
FETI-DP algorithm has been proved independent of the number of subdomains 
for both second-order and fourth-order elliptic systems, and in both two and 
three dimensions, cf. [261 120) . 

The FETI-DP algorithm was first extended to solving incompressible Stokes 
equations by Li [21]. In addition to the coarse level primal velocity variables, the 
subdomain average pressure degrees of freedom are also selected as the coarse 
level variables and the resulting coarse level problem is symmetric indefinite. 
The reduced system for the Lagrange multipliers is still symmetric positive 
definite and a preconditioned conjugate gradient method can be used. Only 
the Dirichlet preconditioner was studied in [2T], and it was proved for two- 
dimensional problems that, under the condition that both the subdomain corner 
and certain edge-average velocity degrees of freedom are selected as coarse level 
primal variables, the condition number bound is independent of the number 
of subdomains and grows only polylogarithmically with the size of individual 
subdomain problems. 

Recently, Kim, Lee, and Park [15] introduced a different FETI-DP formula- 
tion for solving the incompressible Stokes problems, where no pressure variables 
are selected as coarse level primal variables and the resulting coarse level prob- 
lem is symmetric positive definite. Only the lumped preconditioner was studied 
in [IS], for which the edge-average velocity degrees of freedom are no longer 
needed in the coarse level problem; in fact, as few coarse level primal variables 
as for solving positive definite elliptic problems were used and as strong condi- 
tion number bound was established. Reduction of the coarse level problem size 
has also been achieved by Dohrmann and Widlund in an overlapping Schwarz 
type algorithm for solving almost incompressible elasticity, cf. [5,6. Keeping 
the size of the coarse problem small is important in large scale computations; a 
large coarse problem can be a bottleneck and additional efforts in the algorithm 
are needed to reduce its impact, cf. [32 ] 133 ] I3T ] H8 ] 14] 1X6 ] [34] . 

All above mentioned algorithms and their analysis for solving incompress- 
ible Stokes problems are only valid for finite element discretization using dis- 
continuous pressures. Discontinuous pressures have also been used in domain 
decomposition algorithms for other similar type saddle-point problems; see 

[ni m\ us m i m m ee m on . 

In [55], the authors proposed a new FETI-DP algorithm for solving incom- 
pressible Stokes equations, which allows using both discontinuous and contin- 
uous pressures in the finite element discretization. The lumped preconditioner 
was studied in 22. and, similar to [15], as few coarse level primal variables as 
for solving positive definite elliptic problems were used and as strong condition 
number bound was established. 
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The purpose of this paper is two- fold. First, the FETI-DP formulation pro- 
posed in [22] is used as a unified framework to describe the two previous FETI- 
DP algorithms studied in [3T] and [TS]. It is observed that these two FETI-DP 
algorithms can be represented as special cases of using discontinuous pressures 
in the new formulation. The condition number bound estimate based on the 
new formulation also greatly simplifies the analysis in [3T] and [TS]. Second, a 
new Dirichlet preconditioner is studied for the FETI-DP algorithm presented 
under the unified framework using either continuous or discontinuous pressures. 
The same condition number bound as in |21j is obtained. Moreover, this new 
Dirichlet preconditioner involves solving symmetric positive definite subdomain 
problems and is less expensive compared with the Dirichlet preconditioner used 
in [21] where subdomain saddle-point problems need be solved. To stay fo- 
cused on the purpose of this paper, the presentation of the algorithms and their 
analysis is restricted to the case of solving two-dimensional problems. 

The rest of this paper is organized as follows. The finite element discretiza- 
tion of the incompressible Stokes equation is introduced in Section [5] A domain 
decomposition approach is described in Section|3] A reduced system of equations 
is derived in Section [4] Section [5] provides some techniques used in the condi- 
tion number bound estimate. The lumped and the Dirichlet preconditioners are 
studied in Sections [5] and [7] respectively. At the end, in Section [5] numerical re- 
sults for solving a two-dimensional incompressible Stokes problem demonstrate 
the performances of the discussed algorithms and their connections. 



2 Finite element discretization 

We consider the following incompressible Stokes problem on a bounded, two- 
dimensional polygonal domain with a Dirichlet boundary condition, 

!-Au + Vp = f , in fl , 
-V • u = 0, in Q , 
u = uan, on dfl . 

The boundary velocity uga satisfies the compatibility condition J gn ugji ■ n = 0. 
Without loss of generality, we assume that ug^ = 0. 

The weak solution of JTJ is given by: find u E (H^{Q)) 2 = {v e (H 1 ^)) 2 \ 
v = on and p G L 2 (fl), such that, 




o(u,v)+6(v,p) - (f,v), VvG^ffi)) 2 , 
b(u,q) = 0, V<z6L 2 (f>), 



where 

a(u,v)= / Vu-Vv, b(u,q) =- /(V-u)g, (f,v)=/f-v. 

Jn Jn Jn 

The solution of ([2]) is not unique, with p different up to an additive constant. 
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A mixed finite element method is used to solve $Z§. fi is triangulated into 
shape-regular elements of characteristic size h. W 6 (i?g(f2)) represents the 
velocity finite element space and it contains continuous functions. The pres- 
sure finite element space is represented by Q C L 2 (fl). Both continuous and 
discontinuous pressures can be used in our algorithm. A mixed finite element 
space with discontinuous pressure is shown on the left in Figure [1] on a uniform 
triangular mesh of a rectangular domain, where the velocity is piecewise linear 
on the mesh and the pressure is a constant on each union of four triangles as 
shown on the right in the figure. This mixed finite element was used in |21j . 




Figure 1: A mixed finite element with discontinuous pressures. 

A mixed finite element space with continuous pressure is the modified Taylor- 
Hood mixed finite element, as shown in Figure [U cf. [TJ Chapter III, §7]. The 
velocity finite element space contains the piecewise linear functions on the finest 
triangular mesh and the pressure finite element space contains the piecewise 
linear functions on the coarser triangular mesh with the doubled mesh size. 




□ velocity and pressure 
• velocity only 



Figure 2: The modified Taylor- Hood mixed finite clement. 
The finite element solution (u,p) e W0Q of @ satisfies 



(3) 



" A 


B T ' 




u 




f " 


B 







P _ 








Here A, B and f represent respectively the restrictions of a(-, •), &(•, •) and f(-) 
to the finite-dimensional spaces W and Q. In this paper the same notation 
is used to represent both a finite element function and the vector of its nodal 
values. 

The coefficient matrix in ([3]) is rank deficient. A is symmetric positive def- 
inite. The kernel of B T , denoted by Ker(B T ), is the space of all constant 
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pressures in Q. The range of B, denoted by Im(B), is orthogonal to Ker{B T ) 
and is the subspace of Q consisting of all vectors with zero average. The solution 
of @ always exists and is uniquely determined when the pressure is considered 
in the quotient space Q/Ker(B T ). In this paper, when q e Q/Ker(B T ), q 
always has zero average. For a more general right-hand side vector (f , g) given 
in ([3]), the existence of its solution requires that g € Im(B), i.e., g has zero 
average. 

Both mixed finite elements shown on Figures Q] and [2] are inf-sup stable in 
the sense that there exists a positive constant j3, independent of h, such that 

(4) sup f B ™ ] [ > P 2 (q, Zq) , Vq G Q/Ker{B T ). 
wew (w, Aw) 

Here, as always in this paper, (•, •) represents the inner product of two vectors. 
The matrix Z represents the mass matrix defined on the pressure finite element 
space, i.e., for any q G Q, \\q\\^2 — (q,Zq). Z is spectrally equivalent to h 2 I, 
where / represents the identity matrix of the same dimension, i.e., there exist 
positive constants c and C, such that 

(5) ch 2 I <Z< Ch 2 I, 

cf. [351 Lemma B.31]. Here, and in other places of this paper, c and C represent 
generic positive constants which are independent of the mesh size h and H 
(discussed in the following section). 



3 A non-overlapping domain decomposition ap- 
proach 

The domain f2 is decomposed into N nonoverlapping polygonal subdomains fl; , 
i = 1,2, ...,N. Each subdomain is a union of a bounded number of elements, 
with the diameter of each subdomain in the order of H . The nodes on the 
boundaries of neighboring subdomains match across the subdomain interface 
r = (U9fij)\9f2. r is composed of subdomain edges, which are regarded as 
open subsets of T, and of the subdomain vertices, which are end points of edges. 

The velocity and pressure finite element spaces W and Q are decomposed 
into 

W = W ; 0W r , Q = Q/0Qr, 

respectively. Here Wj and Qi are the direct sums of subdomain interior velocity 
spaces Wj , and subdomain interior pressure spaces respectively, i.e., 

N N 

W/ = 0wf, Qz = ©Qi°. 

i=l i=l 

Wr is the subdomain boundary velocity space and it contains the subdomain 
boundary velocity degrees of freedom shared by neighboring subdomains. Qr 
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contains the subdomain boundary pressure degrees of freedom shared by neigh- 
boring subdomains. For the case of using discontinuous pressures, no pressure 
degrees of freedom are shared by neighboring subdomains and Qr is empty. In 
fact, for the discontinuous pressure case, the algorithm presented in this paper 
allows that Qr either be empty or contain any number of intrinsically subdo- 
main interior pressure degrees of freedom; more details on this will be discussed 
in Sections |4~T1 and l4~2l 

To formulate our domain decomposition algorithm, we introduce a partially 
sub-assembled interface velocity space 

Wr = W„ W A = W„ (0 W« \ . 

Here, Wn is the continuous coarse level velocity space and the coarse level 
primal velocity degrees of freedom are shared by neighboring subdomains. The 
complimentary space Wa is the direct sum of subdomain dual interface velocity 
spaces , which correspond to the remaining interface velocity degrees of 
freedom and are spanned by basis functions which vanish at the primal degrees 
of freedom. Thus, an element in the space Wr has a continuous primal velocity 
component and typically a discontinuous dual velocity component. 

In this paper, two choices of Wn are used. In the first, Wn is spanned 
by all the subdomain corner velocity nodal basis functions and the coarse level 
primal variables are only the subdomain corner velocity variables. In the second, 
besides all the subdomain corner velocity nodal basis functions, on each edge 
r u shared by neighboring subdomains fii and Vtj, Wn is also spanned by an 
edge-average finite element basis function such that J rij 1 n y = 0, for any 

£ W^' . Here n,-j denotes a fixed selection of the normal to r y . Therefore 
for the second choice of Wn, the following divergence free boundary condition 

(6) / w«-n = 

is satisfied for all w A <G W A . For more details on choosing coarse level primal 
edge-average velocity variables to satisfy the divergence free condition for in- 
compressible Stokes problems, including for the three-dimensional case, see [Ml 
Section 7] . We note that the choice of Wn depends on the preconditioner used 
in the algorithm. The first choice is sufficient for using the lumped precondi- 
tioner, but for the Dirichlet preconditioner the second one has to be used; for 
more detailed discussions, see Sections [6] and [JJ 

The functions wa in Wa are in general not continuous across T. To enforce 
their continuity, a boolean matrix Ba is constructed from {0, 1, —1}. Each row 
of Ba only contains two non-zero entries, 1 and —1, corresponding to the same 
velocity degree of freedom on each subdomain boundary node, but attributed 
to two neighboring subdomains. For any wa in Wa, each row of Bawa = 
implies that the two degrees of freedom from the two neighboring subdomains 
be the same. When non-redundant continuity constraints are enforced, -Ba has 
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full row rank. We denote the range of Ba applied on Wa by A, the vector 
space of the Lagrange multipliers. 

In order to define a certain subdomain boundary scaling operator, we intro- 
duce a positive scaling factor 8>(x) for each node x on the subdomain bound- 
ary r. Let J\f x be the number of subdomains sharing x, and we simply take 
6'(x) — 1/Af x - In applications, these scaling factors will depend on the heat 
conduction coefficient and the first of the Lame parameters for scalar elliptic 
problems and the equations of linear elasticity, respectively; see [201 [TH] . Given 
such scaling factors on the subdomain boundary nodes, we can define a scaled 
operator Ba,d- We recall that each row of Ba has only two nonzero entries, 1 
and —1, corresponding to the same subdomain boundary node x. Multiplying 
each entry by the scaling factor <$t(x) gives us Ba,d- 

Then solving the original fully assembled linear system ([3]) is equivalent to: 
find ( U/ , Pl , ua, u n , P r, A) e W/ Qj W A W n Qr A, such that 
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Corresponding to the one-dimensional null space of ([3]) , we consider a vector 
of the form (uj, pi, ua, un, pr, A) = (0, l Pl , 0, 0, l Pr ,A), where l Pl and 
l pr represent vectors with value 1 on each entry of the vector. Substituting it 
into (J7J gives zero blocks on the right-hand side, except at the third block 

1, 



(8) 



fA — [Bja B rA ) 



1, 



+ Bl\. 



The first term in ([S]) represents the line integral of the normal components 
of the velocity finite element basis functions across the subdomain boundary, 
and corresponding to the same velocity degree of freedom on the subdomain 
boundary, their values on the two neighboring subdomains are negative of each 
other. Therefore 



[B 



T 
IA 



B r a\ 



1 



pi 



BaBa,d[B 



T 
IA 



B r a] 



1 



pi 



from which we know that fA = in (|SJ), for 

A = — Ba,d[BJ a Bp A ] 
Therefore, a basis of the one-dimensional null space of ([7]) is 



(9) 



0, 1 PI , 0, 0, 1, 



-B a ,d[BJa Bf A ) 



pr 
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which is also valid if pr in (J7J) is empty and then the block l pr in ^ disappears. 

For the case when the coarse level primal velocity space Wn contains both 
the subdomain corner and edge-average variables such that the divergence free 
boundary condition ^ is satisfied, we have Jn . V ■ 



0, for all G w£\ 



which is in matrix form 



iT R (i) 
Pi 



iT R (i) 
Pr 



w«=0. 



As a result 



(10) 



[B. 



T 
/A 



•Br a] 



'-pr 



= 0. 



Therefore, when the divergence free boundary condition ([6]) is enforced and pr 
in ([7J is empty, we know from (I10p and ([9]) that the leading four-by-four diagonal 
block in the coefficient matrix of ([7]) is singular and its null space consists of all 
vectors with a constant pressure and zero velocity. 



4 Reduced system of linear equations 



We first describe the FETI-DP formulation proposed in [35] and then use it as 
a framework to represent the two previous FETI-DP algorithms studied in |21) 
and [TS]. Based on (O, denote 
(11)" 
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The variables (u/, pi, ua, un) can be eliminated from ([7]) and we obtain a 
Schur complement problem for the variables (pr, A) 



(12) 



where 



(13) 



G 



Pr 
X 



G = BcA^Bl, 
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B c A^f. 



Remark 1 The only case that A is singular in the algorithm is when pr in ^ 
is empty and the divergence free boundary condition (0) is enforced, as discussed 
at the end of Section [21 However the definitions in U3\) are still valid, since f 
and columns of B^, are in the range of A and the kernel of A is a subspace of 
the kernel of Be- For the simplicity of notation, we still use A^ 1 in U3\) and in 
other places of this paper to represent the solution of system of linear equations, 
not necessarily the inverse of A. 
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We can see that — G is the Schur complement of the coefficient matrix of (J7J 
with respect to the last two row blocks: 



I 



A 

B c 



B T c 








I 





'a 




-G 



If A is nonsingular, then from the Sylvester law of inertia we can see that 
G is symmetric positive semi-definite and its null space is derived from the null 
space of the original coefficient matrix of 0, cf. , and a basis is given by, 



Ipn -Ba,d[ B Ja B Ta\ 



We denote X = Qr © A. The range space of G, denoted by Rg, is a subspace 
of X. Rg is orthogonal to the null space of G and thus has the form 



(14) R G 



9pr 

9x 



pr ^Pr 



T 

9\ 



Ba,d[bJ a B^ a ] 



When pr in ([7]) is empty and the divergence free boundary condition © is 
enforced, as discussed in Remark [TJ A is singular and then G becomes positive 
definite. The range space formula ([H]) is still valid, cf. pU|) . and in fact Rg 
becomes the whole A. 

In both cases, the restriction of G to its range Rg is positive definite. The 
fact that the solution of ([7]) always exists for any given (f/, fA, fn) on the 
right-hand side implies that the solution of (|12l) exits for any g defined in (fT5|) . 
Therefore g € Rg- When the conjugate gradient method is applied to solve (fT2j) 
with zero initial guess, all the iterates are in the Krylov subspace generated by 
G and g, which is also a subspace of Rg, and where the conjugate gradient 
method cannot break down. After obtaining (pr, A) from solving (fT"2)) , the 
other components (u/, pi, ua, un) in (|7|) are obtained by back substitution. 

The main computation of multiplying G by a vector is the product of A^ 1 
with a vector in the structure of /. We denote 
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and define the Schur complement (coarse level problem) 

Sn = Ann — An r A^A r n- 
Then the product A^ 1 / can be represented by 



(15) 
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which requires solving two subdomain Neumann type problems and one coarse 
level problem. 

The formulation of the reduced system (TT2"j) is valid for using both discon- 
tinuous and continuous pressures in the algorithm, which is discussed in the 
following based on whether pr is empty or not. We note that when the dis- 
continuous pressure is used in the algorithm and pr is empty, (|T2"j) is the same 
system as those obtained in [21] and |15j . 

4.1 pr is non-empty 

This happens when a continuous pressure finite element space is used, where 
Pr represents all the subdomain boundary pressure degrees of freedom shared 
by neighboring subdomains. Then (1121) is a system for both the subdomain 
boundary pressures and the Lagrange multipliers. A in ([T^| and A rr in (TT5|) arc 
both invertible and the coarse level problem operator Sri is symmetric positive 
definite. 

pr can also be non-empty for the case of using discontinuous pressures. Since 
there are no pressure degrees of freedom shared by neighboring subdomains, it 
is free in the algorithm to choose any number of subdomain pressure variables as 
pr- If pr contains at least one pressure degree of freedom from each subdomain, 
then A in (fT"3")) and A rr in (jT5j) are still invertible and Su is symmetric positive 
definite. 

4.2 pr is empty 

pr can be empty only when discontinuous pressures are used in the finite el- 
ement space since no pressure degrees of freedom are shared by neighboring 
subdomains. 

When pr is empty, (|12[) becomes a system for the Lagrange multipliers only 
and it is the same equation as those obtained in [21] and |15j . The implemen- 
tation of the product of A" 1 with a vector as specified in (|15|) is the same as in 
|15j and the resulting Su is symmetric positive definite. Kim et. al. [15] consid- 
ered only the first choice of the coarse level primal velocity space Wn, namely 
it contains only the subdomain corner velocities, for which they proved in [151 
Lemma 3.1] that A and A rr are both invertible. However, A and A rr both be- 
come singular when the second choice of Wn is used to enforce the divergence 
free boundary condition © (required for using the Dirichlet preconditioner, cf. 
Section [7]) , even though their singularities do not affect the multiplication of 
A^ 1 by a vector; see Remark [T] 

In [3T], the divergence free boundary condition ^ is enforced for using 
the Dirichlet preconditioner and a different implementation of multiplying A -1 
with a vector was used to avoid the singularity of A rr in (|15l) . There the 
subdomain constant pressures are pulled out from Qi to form another vector po 
and Qi contains only subdomain interior pressures with zero average on each 
subdomain; po is combined with the coarse level primal velocity variables to 
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form the coarse level problem. More precisely, A in (fTT|) is represented by 
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where, e.g., BJj and B® l7 represent blocks corresponding to the subdomain 
interior pressures with zero average and the subdomain constant pressures, re- 
spectively. Denote correspondingly 
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where A rr is invertible, and dehne the Schur complement (coarse level problem) 
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which is a saddle point problem. Then the product A 1 / can be represented by 
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Remark 2 Even though the implementations of multiplying A" 1 with a vector 
proposed in \21^ and \15^ are different, the same system \l}3jl for the Lagrange 
multipliers is solved. When equipped with the same type preconditioners, their 
convergence rates are the same. 



5 Some techniques 

We first define certain norms for several vector/function spaces. We denote 
(16) W = W 7 0W r . 

For any w in W, denote its restriction to subdomain fii by wW, A subdomain- 
wise -f/^-seminorm can be defined for functions in W by 
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Several of the following lemmas have been proved in [32] ; they are presented 
here for the completeness of this paper. We denote in ([7]) 



(17) 



B = 



Bn Bi& Bin 
Bn BrA -Brn 



The following lemma on the stability of B is [22J Lemma 5.1]. 
Lemma 1 For any w G W and q G Q, (5w, q\ < \~w\ H i WqWl 2 - 



Pmf : (sw, 9 ) 2 =(|/ n v.w«> 5 ) <(|^// a ivw«)p 



We define 



W = w /0Q/0 w A0Wn, 

and its subspace 

(18) W = {w — (w/, p I: wa, w n ) G W : i?//w/ + B/awa + B /n w n = 0} . 

For any w = (w/, pi, wa, wn) G Wo, 
(19) 
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i.e., {-,-)a defines a seminorm on Wq. In (fT9)l . the superscript W is used to 
represent the restrictions of corresponding vectors and matrices to sub-domain 
fii. 

The following lemma is |221 Lemma 6.6]. 
Lemma 2 For any w — (w/, p/, wa, wn) 6 Wo, -Bc^ G -R<3- 

Proof: We know for any (f 7 , f A , f n ) G W, W A W n , <? defined in (p) 
is in For any w = (wj, pj, wa, wn) G Wo, from the definition of A in 
(fTTj) . there always exists (f>, fA, fn) G W/ W A W n , such that 
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Taking such (fj, fA, fn), 9 denned in (fTBf is Bqw. □ 
From (jlip , we denote the first row of Bq by 



Br 



[Br/ i?r a B 



rn 



for the second row, we denote the restriction from W onto Wa by Ra, such 
that for any w = (w/, pj, wa, Wn) € W, ]?aw = wa- Then G, defined in 
(fTU|) , can be written as the following two-by-two block structure 



(20) 



G 



G 



PvPv ^prA 
G\p r G\\ 



where 



Gprpr — By A 



G 



Apr 



B A R A A 



G Pr \ — Br A 1 BaBa, 
G\\ = B A R A A~ 1 R A B A . 



We define a certain jump operator across the sub-domain interface T. Let 
Pd.i ■ W -> W, be defined by 

Bd,i = R^B A d BaRa, 



cf. [23]. We can see that application of Pd,x to a vector essentially computes 
the difference (jump) of the dual velocity components across the sub-domain 
interface and then distributes the jump to neighboring subdomains according 
to the scaling factor 8>{x). In fact, from the definition of Bd,i, the only com- 
ponent involved in its application is the component in the space Wa; all other 
components are kept zero and they are added into the definition to make Pd,i 
more convenient to use in the analysis. 

Note that for any w eW, (Pdaw, Pd.iw) t = (Bl d B a w a , BT d Bawa) . 
The following lemma can be found essentially from Section 6]. 

Lemma 3 There exists a junction <&i(H/h), such that {Pd,iui, Pd.iw) ^ — 
&i(H/h) (w, w)^, for all w £ Wq. For two-dimensional problems, when the 
coarse level primal velocity space contains the subdomain corner velocity vari- 
ables, §i(H/h) = C(H/h)(l + log (H/h)). 

To improve the bound on the jump operator, the jumps across the subdomain 
interface can be extended to the interior of subdomains by subdomain discrete 



harmonic extension. We define a Schur complement : — » by, for 
any w« e w£\ 



(21) 



A d) 

(0 



.4 



A/ 



A (i) 
I A 

(i) 



A 



AA 





r w " 




" 






(*) 

. W A . 









We can see that to multiply H A by the vector w^', subdomain elliptic prob- 



lems with given boundary velocity and w 



,0) 



need be solved. Using 
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H A \ we define the second jump operator Pd.i 



W — » W, by: for any given 
w = (w/, pi, wa, wn) G W, on each subdomain fij, the subdomain interior 



velocity part of Pd,2 w is taken as w^' in the solution of (|2"Tj) . with given subdo- 



main boundary velocity w A = B^'^Bawa. Here B A ' D represents restriction 
of Si. 

D on subdomain fli and is a map from A to . The other components 
of Pd,2W are kept zero. Therefore 



N 



i=\ 
N 

- E 

N 

= E 



(22) 



r (*) " 


T r 


(0 

. W A . 





4« 

(0 

A/ 



.4 



A® 
A IA 

(i) 



A 



AA J L 





r « i 




w} ' 




W 




. W A . 



(*) T _ H -W„,(0 



TV 



a = ^ w A B A B^ D H^B^ D B A w A 

i=i 

2 



B^Bawa 




<$ 2 (H/h)) |w 



1 



where w{?' represents the restriction of (wa, wn) on subdomain f^. The last 
inequality in (|22|) is a well established result, cf., [35j Lemma 6.34]. It has 
been established that, $>2(H/h) = C(l + log (H/h)) 2 , for two-dimensional prob- 
lems, when the coarse level primal velocity space contains the subdomain corner 
velocity variables, cf. [361 Lemma 3.3]. 
Then from {25]) and (JTSJ), we have 

Lemma 4 There exists a function $>2(H/h), such that (Pd,2W, Pd,2w)j < 

<&2{H/h) {w, w}^, for all w G Wo. For two-dimensional problems, when the 
coarse level primal velocity space contains the subdomain corner velocity vari- 
ables, ®2{H/h) = C(l + log (H/h)) 2 . 

The following lemma is also used in our analysis and can be found at [T31 
Lemma 2.3]. 

Lemma 5 Consider the saddle point problem: find (u,p) G W^Q, such that 
(23) 



" A 


B T ' 




u 




' f " 


B 







P 




_ 9 _ 



where A and B are as in j2J), f G W, g G Im(B) C C. 
constant specified in Then \\u\\a < ||f Wa- 1 + \ IMIz- 
matrix defined in Section^ 



Let /3 be the inf-sup 
where Z is the mass 



The lumped and the Dirichlet preconditioners for solving (|12[) will be stud- 
ied in Sections [5] and [7J Let M _1 be symmetric positive definite. When the 
conjugate gradient method is applied to solve the preconditioned system 



(24) 



M~ l Gx = M _1 3, 
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with zero initial guess, all iterates belong to the Krylov subspace generated by 
the operator M~ 1 G and the vector M~ x g, which is a subspace of the range of 
M~ 1 G. We denote the range of M _1 G by Rm~ 1 g- The following two lemmas 
are also from |22j . 

Lemma 6 Let M^ 1 be symmetric positive definite. Then the conjugate gradient 
method applied to solving \2J$ with zero initial guess cannot break down. 



Proof: We just need to show that for any / i £ Rm- 1 Gi Gx ^ 0. Let 
0/i = M^ 1 Gy, for a certain y £ X and y 0. Then Gx = GM^ 1 Gy, which 
cannot be zero since Gy ^ and y T GM~ 1 Gy ^ 0. □ 

Lemma 7 Let M be symmetric positive definite. For any x G R M -i G , 

(Mx,x) — max — — — - -. 

yeRa.y^o (M- 1 y,y) 

Proof: Denote the range of M~^G by R M -i/2 G , then for any A £ Rm- 1 g 

Mix, z ] 

{Mx, x) — ( M5 x, M 5 x ) = max 



□ 



(M?x,M-iy) 2 {y , x y 



yeRa,yjio /M-?y,M-iy) v^cy^o {M l y,y) 



6 The lumped preconditioner 

The lumped preconditioner for solving (fT2|) has been studied in [22] for the 
case of using continuous pressures in the finite element discretization. The 
preconditioner is given by 



Mr 



1 T 



M 



L,\ 



where L pr is the identity matrix of the same length as pr, 



M 



Ba,dRaAR^B^ D 



L,X 



It has been proved in [22.1 that, for the case of using continuous pressures, 
the condition number of the preconditioned operator Mf^G be bounded by 
C$i (H, h)/(3 2 , where &i(H, h) is defined in Lemma [3] and j3 is the inf-sup con- 
stant specified in (jj). The same bound also holds for more general cases where 
pr is non-empty, including for using discontinuous pressures; see Section 14.11 
and [23 Remark 6.10]. 
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When pr is empty, (|T2|) becomes a system for the Lagrange multipliers only, 
which is the same as the one obtained in [15 , cf. Section [4721 The matrix G in 
(fT2)) contains only the second diagonal block in ([20]) . i.e., 

G = G\\ = B A R A A 1 R a B a 

and the lumped preconditioner M^ 1 simplifies to 

Ml 1 = = B A , D R A AR T A B T AD . 

In this section, we establish the condition number bound for the case of 
empty pr following the approach given in [22]. The bound we obtained here 
is the same as that obtained in [15], but this approach greatly simplifies the 
analysis given in [T3]. We also note that this bound is also the same as that 
obtained for non-empty pr in [22] . 

As discussed in [TS] and [22], for using the lumped preconditioner, we only 
need use the first choice of the coarse level primal velocity space, i.e., Wn 
contains only the subdomain corner velocities. We also note that when pr is 
empty, pi contains all the pressure degrees of freedom and the results in Section 
[H] are still valid. In fact the analysis given in the following is essentially the 
special case of the analysis in [22] when the blocks corresponding to pr no 
longer exist. We have the following lemmas. 

Lemma 8 For any w £ Wo, {M^BcW^Bqw) < $ i{H ) h) (^Aw , w^j , where 
&i(H, h) is as defined in Lemma\^ 

Proof: For any given w € Wq, we have 
(M^Bcw^Bcw) = (b a Raw) T M^BaRaw 

= (Ba^aw) T B a ^ d R a AR t a B t ad (Ba^aw) 

(25) = (P D ,iw,P D . 1 w) I <^ 1 {H,h)(w,w) A , 

where we have used Lemma [3] for the last inequality. □ 

Lemma 9 For any given y = g\ 6 Rq, there exits w € Wq, such that Bcw = y, 
and (Aw,w^ < jfr (M^ 1 y 1 y). 

Proof: Given y = g x E R G , let = B T A D g x , u« = (o, u^,o) e 

W,ew A eWn, and = (o, 0, u^, o) 6 W^Q^WAeWn. 

We have 

(26) |u«|^ = (a aa u«u«), 



1G 



and 



(27) 



(i) _ 



[ B A ] 







b a,d9x 




9x, 



where we used the fact that -Ba-BT n = I. 



We consider a solution to the following fully assembled system of linear 

(U) J n ) „(") 



equations of the form ([3]) : find ( u\ , p\ 
that 



uK ; e W/eQ/eWr, such 









A jr " 




r (■"") - 
u / 




" 


(28) 







Bit 




(ii) 
Pi 












Att _ 




(ii) 
u r 








where a particular right-hand side is chosen. We first note that, since g\ G Rq, 
the right-hand side vector of the above system satisfies, cf. (HU), 



It 



-l p ,BiaB a D g\ — 0, 



i.e., it has a zero average, which implies the existence of solutions to the above 
system. Recall that here pi contains all pressure degrees of freedom and there 
are no subdomain boundary pressure degrees of freedom. 

Denote u^' = (vS^ Z \ Up 77 ^ G W. From the inf-sup stability of the origi- 
nal problem ([3]) and Lemma [5l we have 

(29) Iu (//) | 2 ffl <i 



IP 



B/Aui° 



z- 1 



The right-hand side of 



Biau^ 



can be bounded, using Lemma [TJ as follows, 

2 



z- 



= /Bu (7) ,Bu w \ =Cmax 

\ / Z- 1 q&Q 



(30) 



< C max 



|uW|^||g| 

lkll| 2 



i2 =C(A AA u«,u« 



Split the continuous subdomain boundary velocity u r 77 ^ into the dual part 
u^ 7 ' 1 G Wa and the primal part Up 77 ' G Wn, and denote w^- 11 ^ — (uf 1 ^, Pj 77 \ u^ 7 "*, 
We have, from 



(31) 



[ Bu Bja B m ] 



(ii) 
u i 

(ii) 
Pi 

U A 



= -B IA u%> 
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and 



(32) 



B c w (II) = [ B A ] 



(ID 



Pi 



= 0. 



Let w = + wW. We can then see from ([21]) that w e W , cf. (JTHI)- We 
can also see from (|27| and (|32|) that Bc-ii; = Furthermore, by (fl9|) . 



|uW+uW| 2 ffI < luW^ + luW^ < ^(A AA u« ; u«; ; 



where we have used (f2l))) . ([2"5|) and ([3T)|) for the last inequality. 
On the other hand, we have 



A 



□ 



(M^y,y) = gjM^gx = g T x B^ D R A ARlBl D g x = (A AA u^,u 

Theorem 1 For all x = A e R M -i G , 

c/3 2 (M L x,x) < (Gx,x) < $>i(H,h) (M L x,x) , 
where h) is as defined in Lemma\3[ f3 is the inf-sup constant specified in 

Proof: (Gx, x) = x T B c A- 1 B^x = x T BcA^AA^B^x = (Ar x B^x, A^B, 



Since A 1 B]jX G Wq and (-, ■} ? defines an inner product on Wo, we have 



(B c v,x)' 



v, A~ 1 B^x) _ 
(33) {Gx, x) = max ' A = max 

v£W o ,v=£0 \ V 7 V )a oeffc^o (Av,v 



Lower bound: From Lemma [9l we know that for any given y = g\ £ i?c 
there exits w S PFo, such that 5cw = y an d (^Aw,w^ < jj? (Af^ 1 ?/,?/). Then 
from l[3"3")) . we have 



(Gx,x) > 



(Bcw, xY 



> cj3 2 



(y,xY 



' A „..„.) (M L y,y)' 

Since y is arbitrary, using Lemma we have 



(Gx, x) > c/3 2 max 



c/3 2 (M L x, x) 
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Upper bound: From , Lemmas [21 [5] and [71 we have 



(Gx,x) < $i(H,h) max 



(B c v,xY 



< $i(H,h) max 



= $ 1 (H,h)(M L x,x). □ 



v eR G ,y^o (M L L y,y) 

7 The Dirichlet preconditioner 

In the lumped preconditioner discussed in the previous section, the subdomain 
interface jump of the velocity component is extended by zero to the interior of 
subdomains. Comparing Lemmas [3] and SI the discrete subdomain harmonic 
extension of the jump to the interior of subdomains has a better stability, which 
leads to the Dirichlet preconditioner discussed in this section. 

Remark 3 Subdomain discrete Stokes extensions, obtained by solving saddle- 
point problems, are used for the Dirichlet preconditioner studied in \2V§ for solv- 
ing with discontinuous pressures. In the Dirichlet preconditioner proposed 
in this section, the discrete subdomain harmonic extensions of the jump, ob- 
tained by solving symmetric positive definite problems, are used in each iteration. 
Even though these two extensions are spectrally equivalent, using the subdomain 
harmonic extensions is more efficient than solving indefinite subdomain prob- 
lems. 



We define a Schur complement S 



M ( 
any Up = I u 



An 



0w« -> w£'0ws 

is determined by 



W, 



(i) 



(i) 
A - 



ug>)EW«0WS 



r(*) c(») 



^An u r 



For 



(34) 



(i) 
II 

(<) 
II 

(i) 
A/ 

CO 
II/ 

(i) 
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(i) T 
II 



« T 

/A 

« T 
/II 







.4 
B 
.4 
.4 
B 



(i) 
/A 

W 
/A 

(0 
AA 

(i) 
IIA 

(0 

FA 



.4 
B 
.4 
.4 
B 



W 
/ii 

W 
/ll 

(0 
All 

CO 

nn 

(0 
rn 



(i) T 
r/ 



M T 
rA 

(*) T 
i'ii 









r (*) " 
u / 






w 
p} 






(*) 

U A 






(») 
U II 






w 











5 



An u r 

a(i) „(i) 

°An u r 







The solution of (|3"4"f will be needed below in the proof of Lemma [TJ] The 
subdomain incompressible Stokes problem (|34[) contains the subdomain con- 
stant pressure component. To guarantee the existence of its solution, the given 
subdomain boundary velocity Up" 1 = (u^* , Uj^J need satisfy the divergence 

free condition Up • n = 0. In our applications, Up will represent the jump 
of the subdomain interface velocity across the subdomain boundary and the 

(i) 

coarse level component Ujf will always be zero. The divergence free condition 
will be satisfied if ([6]) is enforced. Due to this reason, for the condition number 
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bound analysis for using the Dirichlet preconditioner, we assume that the sec- 
ond choice of the coarse level primal velocity space Wn as described in Section 
[3] is used. 

Remark 4 The assumption that fSJ) is enforced is only required for the condi- 
tion number bound analysis. When (OJ) is not enforced, e.g., when Wn contains 
only the subdomain corner velocity variables, the Dirichlet preconditioner dis- 
cussed below is still symmetric positive definite and the preconditioned conjugate 
gradient method can still be used, even though its condition number bound is no 
longer available. 

The following lemma is a well established result, cf. [2 Theorem 4.1] or [27| 
Lemma 3.1]. 



Lemma 10 For all u^' = [u% , u£> 6 W% Wg ; , 



c #l u r !<?«, - l u rV/»(8(i<) < K l s Wi 
where f3 is the inf-sup constant specified in 0). 

We also denote the direct sum of the discrete subdomain harmonic extension 
operators H$ ,i = l,...,N, defined in (JSTJl, by H A : W A -> W A . 

In the following, the condition number bound for using the Dirichlet pre- 
conditioner is established for the case when pr is non-empty. As discussed in 
Section [4Tj pr in tHj) is non-empty when either the continuous pressure is used 
in the finite clement discretization, or the discontinuous pressure is used and pr 
contains at least one pressure degree of freedom from each subdomain. In fact, 
the same condition number bound also holds for the case when pr is empty, 
which will be discussed briefly at the end of this section. 



We define Mn 1 , by 



(35) M D \ = B A . D H A B AD 

and the Dirichlet preconditioner for solving (TT21 is 



1 T 



Lemma 11 For any w £ W , (M^BcW.Bcw) < C<f> 2 {H,h){Aw,w}, where 
<5>2{H 1 h) is as defined in Lemma^ 

Proof: Given w = (wj, pi, w A , wn) S Wq, let g Pr = BriWi + Br A w A + 
Brnwn- Similar to Lemma El we have from ()22j) . 



(M^Bcw^Bcw) = i (g Pr ,g Pr ) + {B A w A ) T M D \B A w A 
Tp (9pr,9 Pr ) + {B A w A ) T B A:D H A B A ' D B A w A 

(dpr,9pr) + {Pd,2W,P D ,2W)a < ^ (#Pr > 9pr ) + $2(H,h) {w,w) j 
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where we have used Lemma 0] for the last inequality. It is sufficient to bound 
the first term of the right-hand side in the above inequality. 

We denote w = (w/, Wa, Wn) € W. Since B/;w/ + B;AWA + 5/nWn = 0, 
we have 



Biiwi + B/awa + %w n 
-Br/w/ + Bfawa + Brnwn 



_B r /w/ + Bfawa + Brnwn 



(Bw,Bw), 



where -B is defined in (fTT|). From ([5]) and the stability of B, cf. Lemma [TJ wc 
have 



( 3 ^2 (&r , ffpr 



li 2 



Bw,Bw) < C(Bw,Bw 



z- 1 



Cmax 



(Bw,q 



g€Q (q,q). 



q&Q \\q\\ L 2 



where for the last equality, we used the fact that Br/wj + B/a*a + -B/nwn 
and (US). □ 



Lemma 12 Let the coarse level primal velocity space Wn be chosen such that 
(0>P is enforced. For any given y — (g Pr ,g\) 6 Rg, there exits w € Wo, such 

that Bcw — y, and (Aw, w) < jfr (M^y, y). 



Proof: Given y = (g Pr ,g x ) £ R G , let 
each subdomain fli, denote (d/ 



, Pi , Pr 



b a,d9\ and ug 5 = 0. On 
the part obtained through 



the solution of (|3~i)) with given subdomain boundary values = u [ ^' l> and 
uj^ = 0. Let iu^ 1 *) = (uf' l \ p^ % \ u^'^, the corresponding global 

vectors = (u^, p^/\ vS 1 
we know from that w^ € Wo, and 



A ) 



ug 1 , and uO = u£\ uffM. Then 



(37) 



B c w (I) 



Br i 




SrA -Brn 
£?a 



,,(-0 





,9a 



where we used the fact that BaBJ^ d = I. 
Using Lemma [TUl we have 



u 



CM) I? 



\H 1 (Sl i ) 



MA 



< 



C 



U A 
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and summing over the subdomains, 



(38) 



C 



N 



.CM) 
1 A 



We consider a solution to the following fully assembled system of linear equa- 
tions of the form © : hnd(u< 77) , p\ U \ u v U \ p^ 77) ) eW^g^WreOr, 
such that 



(39) 







Air 


B h ' 








- " 


B u 





Bit 







(II) 

Pi 







An 


B Jr 


Att 


B rr 




(ii) 
u r 







. Bpi 





Btt 







C") 
. Pr 




- 9 P T - 



where a particular right-hand side is chosen. Since y G Rq, we can see from 
(fT4")) and (|10[) . that gj r l pr = 0. Therefore the above system has a solution. 

Denote vS 11 ^ — ^u^ 77 \ Up 77 '^. Then from Lemma[S]and ([5]), we have 



(40) 



|„(")|2 , < _L 







< 



z- 1 



C 
fPh? 



Split the continuous subdomain boundary velocity u r 77 ' into the dual part 
and the primal part u[| 7 \ and denote w^- 11 ^ = ^u^ 77 \ p^ 77 \ u^ 7 \ Up 77 
Then we have from that w (7/) G W and 



(41) 



iW 77) 



£?r/ i?rA Brn 
B A 



(II) 

Pi 

„(ii) 



u 



(") 



Let w = w (7) + w( H \ We have w 6 W , B c w = y, from ([37]) and (gT|), and 
from OH) 



C 



I = lu^W 77 '!^ < |u( 7 )|^ + |u( 77 )|^ <^J2 



MA 



c 



H 1 / 2 (dn 



where we have used ([38]) and (|40|) in the last inequality. 
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On the other hand, we have from ([33]) 



(M^y, y) = (g Pr ,g Pr }+ gjM D \g x 



N 

^2 (dvr'dpr) + y ] 



1 



fJX 



CM) 

A 



□ 



H 1 / 2 (dfl i 



With Lemmas QT] and [121 similar to the proof of Theorem [TJ we have the 
following theorem. 

Theorem 2 Let the coarse level primal velocity space Wn be chosen such that 
0) is enforced. For all x — (pr. A) E R M -i G , 

C0 1 (M D x, x) < (Gx, x) < C$ 2 (H, h) (M D x, x) , 

where $2 (H, h) is as defined in Lemma [7J /3 is the inf-sup constant specified in 
&■ 

When pr in (|12p is empty, the matrix G contains only the second diagonal 
block in (|20j) and the Dirichlet preconditioner for (IT21) becomes only » , as 
given in (|35|) . i.e., 

M D X = Md.a = B AiD H A Bl iD . 



Counterparts of Lemmas [TT1 and [T2l for this case can be proved as well; in fact, 
their proofs are essentially the first half in the proofs of Lemmas [TT] and fl2l 
respectively. The condition number bound in Theorem[2]can then be established 
for this case in the same way. 



8 Numerical experiments 



We consider solving the incompressible Stokes problem (UJ in the domain f2 = 
[0, 1] x [0, 1]. Zero Dirichlet boundary condition is used. The right-hand side 
function f is chosen such that the exact solution is 



sin 3 (7rx) sin 2 (7rj/) cos(7ry) 
— sin 2 (7ra) sin 3 (7ry) cos(7r;c) 



and p 



Two mixed finite element discretizations, as shown on Figures Q] and [2] in 
Section [21 are used for the cases of discontinuous and continuous pressures, 
respectively. The preconditioned system ([24]) is solved by the preconditioned 
conjugate gradient method; the iteration is stopped when the L 2 — norm of the 
residual is reduced by a factor of 10~ 6 . 

In each of the following tables, we present the performance of three differ- 
ent variants of the FETI-DP algorithm represented under the same framework, 
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as discussed in Sections 14.11 and 14.21 "continuous pressure" for the case when 
the continuous pressure is used in the algorithm and pr contains all the sub- 
domain boundary pressure degrees of freedom; "discontinuous pressure" for the 
case when the discontinuous pressure is used and pr contains just one pressure 
degree of freedom from each subdomain; "pr empty" for the case when the dis- 
continuous pressure is used and pr is chosen empty. For each case, the extreme 
eigenvalues and the iteration count for each experiment are shown. The two 
methods discussed in [2T] and [TS] solve the same system (TT2")) . and their only 
difference is in the implementation of multiplying G by a vector, cf. Remark [3J 
Therefore their convergence rates are the same, when they are equipped with 
the same type preconditioner, and their performance is reported under the case 
u pr empty" in the tables. 

Tables [T] and [2] show the performance of using the lumped preconditioner for 
different cases with two choices of the coarse level prime variables. In Table [I] 
the first choice is used, namely only the subdomain corner velocities are taken as 
the coarse level primal variables. We can see that for each variant of the FETI- 
DP algorithm, the convergence rate is independent of the number of subdomains 
for fixed H/h; for fixed number of subdomains, the condition number grows 
presumably in the order of (H/h) (1 + \og(H/h)) as established in Section [51 
In Table [2J we test the lumped preconditioner with the second choice of the 
coarse level primal space, namely it contains both the subdomain corner velocity 
variables and the edge-average velocity variables such that ([5]) is enforced, as 
discussed in Section [3| Even though the edge-average velocity variables are not 
required in the coarse level primal space for the lumped preconditioner case, 
including them improves the convergence rate for each method. We also observe 
from Tables Q] and [2] that performances of the three variants of the FETI-DP 
algorithm are quite similar, while the convergence when using discontinuous 
pressure and choosing pr empty is a little faster than the other two cases. 

Tables [3] and [4] show the performance of using the Dirichlet preconditioner. 
In Table [31 only the subdomain corner velocities are taken as the coarse level 
primal variables, for which the divergence free boundary condition ([5]) is not 
satisfied and no scalable condition number bound of the FETI-DP algorithm is 
available. Indeed Table[3]shows that for each variant of the FETI-DP algorithm, 
the convergence rate deteriorates with the increase of the number of subdomains 
when the subdomain problem size is fixed. In Table [4l both the subdomain 
corner and the edge-average velocity variables are taken as the coarse level 
primal variables such that ([6]) is enforced. We can see that for each variant 
of the FETI-DP algorithm, the convergence rate is independent of the number 
of subdomains for fixed H/h; for fixed number of subdomains, the condition 
number grows presumably in the order of (1 + \og(H/h)) as established in 
Section [7] 
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